Identification of miRNA Associated with Linolenic Acid Content Traits in Xiangqingcai

Wang Zhen ^1*　Liu Zhaokun ²　Du Yicheng ¹　Chen Sujuan ¹　Chen Guoyuan ¹

1 Suzhou Polytechnic Institute of Agriculture, Suzhou, 215008; 2 Suzhou Academy of Agricultural Sciences, Suzhou, 215008

* Corresponding author, 18962103126@189.com

Abstract　The purpose of this study is to explore the miRNA and its target genes that affect linolenic acid content in Xiangqingcai (Brassica campestris ssp.chinensis L. Makino, NHCC), and to clarify its molecular function in the regulation of this character. In this research, using the mature leaves of ‘xiuhuajin’ and ‘Heiyexiangqingcai’ which have largest difference in linolenic acid content as materials, the differentially expressed miRNA was identified by high-throughput microRNA sequencing, and the target gene function was predicted and further function annotated. In total 236 miRNAs were detected, of which 159 were known miRNA and 77 were novel ones. In these miRNAs, 28 were up-regulated and 35 were down regulated. A total of 6041 target genes corresponding to the differentially expressed miRNA were annotated, including many fatty acid metabolism pathways, such as "α-linolenic acid metabolism", "fatty acid synthesis", "unsaturated fatty acid biosynthesis" and "fatty acid degradation", etc. The target genes of differential expressed miRNA mainly are some key enzymes of linolenic acid synthesis or metabolism, such as lipoxygenase (LOX), 3-ketoacyl-CoA thiolase (KCS), which play an important role in the synthesis and transformation of linolenic acid. Through negatively regulating their expression and other key genes of fatty acid metabolism and the related genes of linolenic acid synthesis such as KCS through the differential expression of these miRNAs in different varieties, on the one hand reduces the metabolism loss of linolenic acid, on the other hand increases the synthesis amount of linolenic acid, thus resulting in the increase of linolenic acid content.

Keywords　Xiangqingcai, Linolenic acid, Differentially expressed miRNA, MiRNA target genes

香青菜是苏州地区地方传统特色蔬菜，栽培历史悠久，是一种具有浓郁油香的不结球白菜(Brassica campestris ssp.chinensis L. Makino, NHCC) (陈虎根等, 2004)。其下有4个亚种类型，分别是‘黑叶香青菜’、‘黄叶香青菜’、‘绣花筋’和‘黑杂一号’香青菜等，‘绣花筋’香味较淡，而另三者香味浓郁(杨雪梅等，2007)。目前为止，这种奇特的香味只在香青菜中发现含有，然而，有关香青菜的各类研究报道还很少，大多停留在应用传统育种方法进行品种复壮、选育等方面，对于香青菜香味特别是亚麻酸含量的分子机理研究还处于起步阶段(杨雪梅等, 2012; 陈国元和张磊, 2016)。王桃云等(2017)对香青菜香味来源的挥发性芳香物质进行提取，并分析其成分。分析显示，香青菜香味主要化学成分为脂肪族、芳香族和萜类3大类化合物及其含氧衍生物。本人对香青菜通过GC-MS分析发现，‘黑叶香青菜’亚麻酸含量低，而‘绣花筋’亚麻酸含量高达近50%。迄今为止，还未见关于香青菜亚麻酸含量相关的基因功能报道。

研究miRNA在植物的生长发育以及胁迫应答机制中的作用已逐渐成为热点(Cuperus et al., 2011)。miRNA是广泛岑在于生物体的普遍存在的非编码小分子RNA，其长度约20~24 nt。在植物中，miRNA对发育进程、抗逆反应等起着重要的作用(Cho and Paszkowski, 2017; Luo et al., 2018)。在亚麻酸的生物合成中，miR158、miR159、miR166、miR167、miR172调控脂肪酸去饱和酶的表达(王桃云等, 2017)，而miRNA156、miRNA172、miRNA393等通过调控LOX基因的表达，降低亚麻酸的含量(Qiu et al., 2015)。

本研究对‘黑叶香青菜’和‘绣花筋’的成熟叶片进行高通量miRNA测序，鉴定两者表达量存在差异的miRNA并通过软件预测对应的靶基因，筛选出与亚麻酸代谢相关基因的同源基因，如脂氧合酶(LOX)、转氨酶和丙酮酸脱氨酶等(Hu et al., 2014)，从而得到参与香青菜亚麻酸含量调控的miRNA及其对应靶基因，揭示香青菜亚麻酸含量甚至香味差异的分子机制。

1结果与分析

1.1测序数据分析

对测序技术得到的数据信息进行质量管理控制，全部6个样品的clean Data都超过11.99 M，碱基数均大于0.6 G (表1)。经过去除带接头和低质量的序列，最终得到的3个‘绣花筋’样品的Clean Reads分别为11.73 M、14.97 M和12.97 M，而‘黑叶香青菜’的Clean Reads则分别为13.94 M、12.87 M、12.54 M。两个品种的Q30 (质量值大于30)值最低为97.26%和96.90%。

1.2 miRNA的鉴定

在miRBase (V21)中，将以上定位到参考序列上的reads，和十字花科芸苔属序列进行比对，鉴定到的已知miRNA指的是Reads与库中已知miRNA完全相同的序列，得到‘绣花筋’和‘黑叶香青菜’库序列匹配上miRNA的详细情况，其中包括miRNA的序列信息、长度、二级结构等。结果共得到miRNA前体共147，成熟的miRNA159个(表2)。具有发夹结构是miRNA前体的特殊标志，miREvo (Wen et al., 2012)和mirdeep2 (Friedländer et al., 2012)等预测软件可鉴于此预测到新的miRNA，对于没有比对到已知miRNA的sRNA序列，我们利用这两个软件来预测新miRNA，结果共鉴定出新的miRNA前体98个，成熟miRNA共77个(成熟体命名格式按照Bra-miR-n*来, 其中*为数字编号)。图1是部分鉴定得到的已知miRNA和新miRNA的二级结构示意图。

1.3 miRNA家族分析

高通量测序得到的159个已知miRNA分别聚集到30余个miRNA基因家族中，每个家族拥有的miRNA的数目和表达量各有不同。图2显示了本试验鉴定出的有成员超过3个的miRNA家族，miR159、miR164等4个家族在香青菜中都鉴定出3个成员，miR157、miR167等6个家族都鉴定出4个成员，miR169有6个成员，而miR156最多，有7个成员被鉴定出(图2)。其他miRNA家族鉴定出的已知miRNA基本为1~2个，可能具有成员更多的miRNA家族在香青菜的香味调控中发挥更大的作用。

1.4差异表达miRNA的鉴定分析

以miRNA表达分析中的readcount作为源数据，随后使用DESeq2进行基于负二项分布的方法分析来鉴定两个样品间差异表达的miRNA (Love et al., 2014)。试验共筛选出两个样本的差异表达miRNA共63个，其中，以‘黑叶香青菜’作为对照，在‘绣花筋’中上调的miRNA有28个，下调的有35个。已知的miRNA中，有35个miRNA在两个品种中出现差异表达，其中上调表达的有12个，下调表达的有23个；新的miRNA中，共28个显示差异表达，其中上调表达的miRNA有16个，下调表达的miRNA有12个(表3)。

1.5靶基因预测

通过鉴定到的miRNA分别与其对应植物物种中基因的碱基数据，使用psRNATarget进行靶基因的预测，预测到靶基因的miRNA共有234个，靶基因数量为17 036个。这些miRNA中，有159个是已知的miRNA，预测到靶基因13 128个；有75个是新的miRNA，预测到靶基因有3 908个(表4)。

1.6香青菜差异表达miRNA的靶基因功能注释

1.6.1 GO注释

对上述得到的靶基因进行功能比对注释，在注释库中得到了63个miRNA对应的注释，共计靶基因6 041个，包括GO注释库中4556个差异表达miRNA的靶基因和KEGG注释库中的2 525个靶基因。

通过进行GO在线分析，研究一共获得了4 556个差异表达miRNA的靶基因被注释到。这些靶基因共定位到GO的29个功能群当中，分别隶属于两个大的功能亚类，包括生物学过程亚类(biological process，BP)9个和分子功能亚类(molecular function, MF)20个，分析中并没有发现与细胞组份相关的差异miRNA的靶基因。在生物学过程中，“蛋白质磷酸化”占比超过40%，为最大的功能群，其次是“转录水平调控”、“酰胺代谢相关调控”、“转录后调控”、“蛋白质以及代谢过程调控”、“细胞蛋白代谢调控”等，占比均在10%左右，而“病毒诱导的寄主细胞融合”则占比最低，仅占1%左右(图3)。分子功能分类是差异表达miRNA靶基因注释到的最多的一类，除“转录调控活性”、“ADP结合”和“ATP合酶活性”这三类外，其他17类得到注释的靶基因均超过800个，而“水解酶活性”占比最多。

1.6.2 KEGG通路分析

用KEGG对差异表达的miRNA的靶基因进行分析，发现其参与的主要通路有20条，其中，“植物病原相互作用”通路有49个靶基因注释到，占比11.3%，41条靶基因注释到“自噬调节”通路，占比9.3%，31条靶基因注释到脂肪酸代谢中，占比7.0%，“过氧化物酶体”通路有28条靶基因注释到，其余包括“氰氨酸代谢”、“磷酸肌醇代谢”、“核苷酸切除修复”等通路都有20条以上靶基因得到注释。值得注意的是，一些和脂肪酸代谢相关的通路，如“乙醛酸和二元酸代谢”、“脂肪酸降解”、“脂肪酸合成”、“不饱和脂肪酸的生物合成”、“脂肪酸延长”等都有10条以上基因得到注释，而与亚麻酸含量直接相关的“α-亚麻酸代谢”通路也有15条靶基因注释(图4)。

1.7香青菜中与亚麻酸合成相关的miRNA

根据靶基因功能注释的结果，对功能为“α-亚麻酸代谢”的靶基因及其对应的差异表达miRNA进行重点分析，共发现15个靶基因及其对应的15个差异表达miRNA，包括bra-miR5711、bra-miR-n62、bra-miR9555a-3p、bra-miR-n165、bra-miR-n34、bra-miR158-5p、bra-miR397a、bra-miR845a、bra-miR390a-5p、bra-miR156a-3p、bra-miR5716、bra-miR824-5p、bra-miR-n136、bra-miR161-5p和bra-miR860-3p等。其中，将‘黑叶香青菜’作为对照时，‘绣花筋’中6个miRNA的表达上调，上调幅度最大的是bra-miR-n136，对应靶基因为LOX6同源基因；下调的有9个miRNA，下调幅度最大的是bra-miR5716，对应靶基因为ACX-4同源基因。

2讨论

植物的脂肪酸代谢系统中，脂氧合酶(LOX)将亚油酸和亚麻酸是作为主要底物，两者被其催化反应生成氢过氧化脂肪酸，经过氢过氧化物裂解酶(HPL)的裂解作用继续转化为C6、C9醛类，随后通过乙醇脱氢酶(ADH)的氧化反应形成对应的醇类物质，酰基转移酶(AAT)参与醇类向酯类的转化(Inés et al., 2015; Kim et al., 2015)。因而LOX基因的表达水平直接影响亚麻酸的含量。在本研究中，有4个LOX基因(LOX2、LOX3、LOX4和LOX6)作为差异表达miRNA的靶基因影响香青菜的亚麻酸含量，其对应miRNA为bra-miR5711、bra-miR-n62、bra-miR390a-5p、bra-miR156a-3p、bra-miR-n136等。在这5个miRNA中，虽然bra-miR-n62和bra-miR156a-3p在‘绣花筋’中表达量下调(log2倍数分别是-1.03和-0.76)，但其他3个均为上调且幅度较大，bra-miR-n136更是只在‘绣花筋’中鉴定出，为其特有的miRNA，其readcounts高达172.23，相较于‘黑叶香青菜’，可靶向LOX6基因并大幅度降低其表达量，再结合bra-miR390a-5p的上调对于LOX2的负调控，表明了这种miRNA的调控机制可能对香青菜亚麻酸的含量多寡有十分关键的作用。

此外，bra-miR9555a-3p和bra-miR-n-165的靶基因功能注释为酰基辅酶A氧化酶(Acyl-CoA oxidase，ACX)，其参与脂肪酸代谢下游醇类到酯类的转化(Adham et al., 2005)，后者是桃果实香味产生的主要来源，ACX是该途径的关键限速酶(Arent et al., 2008; 齐玉洁等, 2012; Fan et al., 2017)。bra-miR9555a-3p和bra-miR-n-165在‘绣花筋’中的上调表达将抑制ACX基因的表达，引起代谢过程中间产物的积累，最终氢过氧化脂肪酸的累积将抑制LOX催化的亚麻酸的转化，从而增加其在‘绣花筋’中的含量。脂肪酸代谢中，上游关键基因LOX和下游关键基因ACX的表达水平降低，使得亚麻酸更少的被代谢消耗，而醇类、酯类等芳香物质的较少在表观上导致了‘绣花筋’的香味较‘黑叶香青菜’更淡。

差异表达miRNA中，bra-miR845a靶向4个基因，而后者都为3-酮酰基-CoA硫酶2基因(3-ketoacyl-CoA thiolase, KCS)同源基因。KCS是超长链单不饱和脂肪酸(VLCMFA)生物合成途径的关键酶，而VLCMFA的合成前体主要以油酸为主(Deyu et al., 2015)。Liu和Li的(2014)报道表明，油酸是亚麻酸的重要合成前体。因而，VLCMFA和亚麻酸存在对油酸这一共同合成前体的竞争关系，当VLCMFA合成增多时，必然无法产生较多的亚麻酸。在‘绣花筋’小RNA文库中，bra-miR845a的readcounts为所有鉴定出的差异表达miRNA中最高的，达到1 223.13，而‘黑叶香青菜’中该值为286.99。这一miRNA的大幅表达水平差异将导致在‘绣花筋’中KCS基因的表达大幅抑制，从而减少了油酸向VLCMFA的生物合成，使得亚麻酸的合成有充足的前体使用。

本研究中，通过对亚麻酸含量差异较大的两个香青菜品种的高通量测序分析，得到15个和该性状高度相关的miRNA及其对应的15个参与“α-亚麻酸代谢”通路的靶基因，其中有3个miRNA是新发现的(bra-miR-n62, bra-miR-n165和bra-miR-n136)。推测香青菜通过对这些miRNA在不同品种中的差异表达负调控LOX、ACX等脂肪酸代谢关键基因和KCS等亚麻酸合成相关基因的表达，从而一方面减少亚麻酸的代谢损失，另一方面增加亚麻酸的合成量，最终导致亚麻酸含量的增高(‘绣花筋’中)。本试验初步揭示了miRNA及其靶基因参与的香青菜亚麻酸含量性状调控的分子机理，为更好的研究和利用高含量亚麻酸品种的香青菜，更好的开发其保健作用打下理论基础。

3材料与方法

3.1试验材料

本试验所选用的材料为苏州蔬菜研究所提供的多年自交‘绣花筋’和‘黑叶香青菜’纯合种。2018年9月种植于苏州农业职业技术学院实训田，分别取其成熟期叶片提取总RNA，每个品种取3组样品设为生物学重复。

3.2试验方法

3.2.1总RNA的提取及smallRNA文库构建

RNA提取使用TRIzlo(Invitrogen)试剂盒，使用15%的TBE-Urea胶分离得到小RNA。纯化的小RNA 连上5’和3’接头(大小为70~90 nt)，设计接头引物并用其进行反转录PCR计15个循环，纯化后的扩增片段进行上机测序(Illumina/Solexa G1 sequencer)。使用Small RNA Sample Pre Kit试剂盒对合格的样品构建cDNA文库，利用小RNA的5’端和3’端分别有完整的磷酸分子基团和羟基的特征，以总RNA为初始材料，将Small RNA两端同时加上5’和3’接头(大小为70~90 nt左右)，反转录PCR15个循环反应后合成cDNA。再经过PCR进行分析扩增，PAGE胶电泳方法分离目标的DNA片段，最后切胶回收得到两个不同材料共6个样品cDNA文库(‘绣花筋’: XHJ1-3, ‘黑叶香青菜’: HY1-3)。

3.2.2 miRNA的鉴定与分析

原始reads经过去接头和过滤低质量后，将cluster reads依次与Rfam库、miRBase库、RepBase库等进行比对，获得小RNA的注释信息，并得到已有的，以sRNA注释信息此作为预测新的miRNA的基础(Griffiths-Jones et al., 2008)。已知的miRNA鉴定分为两步：首先未注释的sRNA序列先与白菜miRNA进行比对，鉴定是否有已知的miRNA，此过程不允许碱基错配；随后将未配对的剩余序列与植物所有miRNA进行比对，比对在其他植物中是否有已知miRNA。若sRNA没有得到任何比对结果，则再通过mireap软件来预测新的miRNA。检测sNRA上下游100 bp左右序列的Dicer酶切位点信息和能量特征，miRNA的二级结构通过RNAfold工具预测得到。

3.2.3差异表达miRNA的筛选及其靶基因功能注释

差异表达miRNA的筛选通过IDEG6软件进行分析，依据fold-change(log2)>1或者fold-change(log2)<-1，且P-value<0.01鉴别miRNA在两个文库中是否存在显著表达差异。第一步是对两个样本中miRNA的表达量进行标准化处理(tags per million, TPM)，然后用TPM计算fold-change和P-value。利用在线软件psRNATarget对鉴定出的miRNA进行预测，得到其对应的靶基因信息。随后对差异表达miRNA的靶基因进行在线功能注释，使用的数据库为主要为GO (Gene Ontology, http://www. http://geneontology.org/)和KEGG (Kyoto Encyclopedia of Genes and Genomes, https://www.kegg.jp/) (郭强等, 2019)。

作者贡献

王镇和刘照坤是本研究的试验设计者和执行人；杜奕承是参与数据整理和论文初稿的写作；杜奕承、陈素娟、陈国元参与部分试验；王镇是项目的负责人，指导试验设计，数据统计，论文写作与修改。

致谢

本研究由江苏省高等学校自然科学研究面上项目(NO: 18KJB210010)和苏州市科技计划农业科技创新项目(NO: SNG2017057)共同资助。

参考文献

Adham A.R., Zolman B.K., Millius A., and Bartel B., 2005, Mutations in arabidopsis acyl-coa oxidase genes reveal distinct and overlapping roles in 2-oxidation, Plant J., 41(6): 859-874

Arent S., Pye V.E., and Henriksen A., 2008, Structure and function of plant acyl-coa oxidases, Plant Physiol. Biochem., 46(3): 292-301

Chen H.G., Yang X.M., Jiang S.D., and Wang X.Y., 2004, The local characteristic variety of Suzhou -Xiangqingcai, Shanghai Shucai (Shanghai Vegetables), (5): 22-23 (陈虎根，杨雪梅，蒋树德，王香云，2004，苏州地方特色品种-香青菜, 上海蔬菜, (5): 22-23)

Cuperus J.T., Fahlgren N., and Carrington J.C., 2011, Evolution and Functional Diversification of MIRNA Genes, Plant Cell, 23(2): 431-442

Cho J., and Paszkowski J., 2017, Regulation of rice root development by a retrotransposon acting as a microrna sponge, eLife, 6: e30038

Chen G.Y., and Zhang L.2016, Purification and rejuvenation of a new variety of Suzhou Dongshan Xiangqingcai, Changjiang Shucai (Changjiang Vegetables), (13): 21-22 (陈国元, 张磊, 2016, 提纯复壮新品种苏州东山香青菜, 长江蔬菜, (13): 21-22)

Deyu T., Shian W., Lihao W., Jialin W., and Fuli L., 2015, The biosynthesis and metabolic engineering of very long-chain monounsaturated fatty acid, Biotechnol. Bulletin, 12: 80-84

Fan J., Yu L., and Xu C., 2017, A central role for triacylglycerol in membrane lipid breakdown, fatty acid β-oxidation, and plant survival under extended darkness, Plant Physiol., 174(3): 1517

Friedländer M.R., Mackowiak S.D., Na L., Wei C., and Nikolaus R., 2012, Mirdeep2 accurately identifies known and hundreds of novel microrna genes in seven animal clades, Nucleic Acids Res., 40(1): 37-52

Griffiths-Jones S., Saini H.K., Van D.S., and Enright A.J., 2008, Mirbase: tools for microrna genomics, Nucleic Acids Res., 36(Database issue): D154

Guo Q., Wang Y.Z., Chen J.J., Zhou L., and Xu B., 2019, Identification of microRNA related to cytoplasmic male sterile line of Alfalfa, Caodi Xuebao (Journal of Grassland), 27(5): 1138-1146 (郭强, 王英哲, 陈晶晶, 周仂, 徐博, 2019, 紫花苜蓿细胞质雄性不育系相关microrna的鉴定, 草地学报, 27(5): 1138-1146)

Hu T., Zeng H., Hu Z., Qv X., and Chen G., 2014, Simultaneous silencing of five lox genes increases the contents of alpha linolenic acid and linoleic acid in tomato (solanum lycopersicum L) fruits, J. Agric. Food Chem., 62(49): 11988-11993

Inés Ponce de León, Hamberg M., and Castresana C., 2015, Oxylipins in moss development and defense, Front. Plant Sci., 6: 483

Kim S.U., Kim K.R., Kim J.W., Kim S., and Park J.B., 2015, Microbial synthesis of plant oxylipins from γ-linolenic acid through designed biotransformation pathways, J. Agric. Food Chem., 63(10): 2773-2781

Liu L.Z., and Li J.N., 2014, QTL Mapping of Oleic Acid, Linolenic Acid and Erucic Acid Content in Brassica napus by Using the High Density SNP Genetic Map, Scientia Agricultura Sinica, 47(1): 24-32

Love M.I., Huber W., and Anders S., 2014, Moderated estimation of fold change and dispersion for rna-seq data with deseq2, Genome Biol., 15(12): 550

Luo M., Gao Z., Li, H., Li, Q., Zhang, C., and Xu, W., 2018, Selection of reference genes for mirna qrt-pcr under abiotic stress in grapevine, Sci. Rep., 8(1): 4444

Qi Y.J, Zhang X., Yang X., Gao Z.S., and Jia H.J., 2012, Cloning and expression analysis of acyl coenzyme A oxidase gene of "hujingmeilu" peach fruit, Zhongguo    Nongye Kexue (Agricultural Science in China), 45(9): 1758-1765 (齐玉洁, 张鑫, 杨夏, 高中山, 贾惠娟, 2012, ‘湖景蜜露’桃果实酰基辅酶a氧化酶编码基因的克隆及表达分析, 中国农业科学, 45(9): 1758-1765)

Qiu Z., Li X., Zhao Y., Zhang M., Wan Y., Cao D., Lu S., and Lin J., 2015, Genome-wide analysis reveals dynamic changes in expression of micrornas during vascular cambium development in chinese fir, cunninghamia lanceolate, J. Exp. Botany, 66(11): 3041-3054

Wang T.Y., Jiang W.N., Hu C.Y., Li Q., Shen X.L., and Lv H.C., 2017, Study on the extraction, chemical composition and antibacterial activity of the essential oil from Xiangqingcai, Zhongguo Liangyou Xuebao (Chinese Journal of Grain and Oil), 32(3): 81-87 (王桃云, 蒋伟娜, 胡翠英, 李倩, 沈雪林, 吕海超, 2017, 香青菜挥发油提取及化学成分和抑菌活性研究, 中国粮油学报, 32(3): 81-87)

Wen M., Shen Y., Shi S., and Tang T., 2012, Mirevo: an integrative microrna evolutionary analysis platform for next-generation sequencing experiments, BMC Bioinformatics, 13(1): 140

Yang X.M., Jiang S.D., Yin Y.L., Dong J.M., and Xu M., 2007, Breeding of common Chinese cabbage Heiza-1. Zhongguo Shucai (Chinese Vegetables), 1(4): 29-30 (杨雪梅, 蒋树德, 尹渝来, 董建明, 徐溟, 2007, 普通白菜(香青菜)黑杂-1号的选育, 中国蔬菜, 1(4): 29-30)

Yang X.M., Han J.J., Zhang S., and Chen H.G., 2012, Breeding of new cabbage variety hanxiangqingcai, Henan Nongye Kexue (Henan Agricultural Sciences), 41(11): 114-116 (杨雪梅, 韩建军, 张胜, 陈虎根, 2012, 小白菜新品种寒香青菜的选育, 河南农业科学, 41(11): 114-116)